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Abstract. This paper describes the recently developed mixed mimetic spectral element method 
for the Stokes problem in the vorticity-velocity-pressure formulation. This compatible discretiza- 
tion method relies on the construction of a conforming discrete Hodge decomposition, that is 
based on a bounded projection operator that commutes with the exterior derivative. The projec- 
tion operator is the composition of a reduction and a reconstruction step. The reconstruction in 
terms of mimetic spectral element basis-functions are tensor-based constructions and therefore 
hold for curvilinear quadrilateral and hcxahcdral meshes. 

For compatible discretization methods that contain a conforming discrete Hodge decompo- 
sition, we derive optimal a priori error estimates which are valid for all admissible boundary 
conditions on both Cartesian and curvilinear meshes. These theoretical results are confirmed 
by numerical experiments. These clearly show that the mimetic spectral elements outperform 
the commonly used /f(div)-compatiblc Raviart-Thomas elements. 



1. Introduction 

Let J] C R", n > 2, be a bounded contractible domain with boundary T = dtt. On this domain 
we consider the Stokes problem, consisting of the equations for conservation of momentum and 
for conservation of mass, 

(1.1a) V • a = f on Q, 

(1.1b) divu = g on f2, 

where the stress tensor a is given by 

(1.2) a = -uVu+pI, 

with u the velocity vector, p the pressure, / the forcing term, g the mass source and v the kinematic 
viscosity. For analysis purposes we choose v = 1. 

This paper considers the recently developed mixed mimetic spectral element method (MMSEM) 
[40, 41]. This compatible finite/spectral element method is based on the compatible discretization 
of the exterior derivative d from differential geometry, which represents the vector operators, 
grad, curl and div. The Stokes problem expressed in terms of these vector operations is known as 
the vorticity-velocity-pressure (VVP) formulation, [9, 23]. For the VVP formulation, the Laplace 
operator is split using the vector identity, —Au — curl curl* u — grad* divu, and by introducing 
vorticity as auxiliary variable, co = curl* u. The VVP formulation of the Stokes problem becomes 

(1.3a) uj — curl* u = 0, on fi 

(1.3b) curio; — grad* divu + grad* p — f, on il 

(1.3c) divu — g, on fi. 
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Following [9, 40] we make a distinction between the operators grad, curl and div, that correspond to 
the classical Newton-Leibnitz, Stokes circulation and Gauss divergence theorems, and the operators 
-grad*, curl* and -div* that are their formal Hilbert adjoints, 



The distinction between the two types of differential operators is made explicitly, because the 
construction of our conforming finite element spaces relies on the three mentioned integration 
theorems, while the mixed formulation relies on the formal Hilbert adjoint relations. While in 
vector calculus this distinction is not common, in differential geometry these structures naturally 
appear since they make a clear distinction between metric-free (topological) and metric-dependent 
operations. 

The MMSEM is a compatible discretization method that relies on the construction of a con- 
forming discrete Hodge-decomposition, which implies a discrete Poincare inequality. It requires 
the development of a bounded projection operator that commutes with the exterior derivative. 
The bounded projection is a composition of a reduction by means of integration and mimetic 
spectral element basisfunctions as reconstruction. 

The reduction onto /c-dimensional submanifolds result in the discrete unknowns representing 
integral quantities. This is one of the major differences with related methods as the Marker and 
Cell scheme [32] and the lowest-order Raviart-Thomas and Nedelec compatible finite elements 
[45, 48], where use is made of averaged quantities. 

The basis functions, used for the reconstruction, are constructed using tensor products of one 
dimensional nodal and edge interpolation basis functions [28] , and therefore hold for quadrilateral 
and hexahedral meshes. They belong to the class of compatible finite elements, and were con- 
structed based on the mimetic framework first described in [38] and later extended in [11]. The 
mimetic framework, including the mimetic spectral elements, were extensively described in [41]. 
This mimetic framework relies on the languages of differential geometry instead of vector calculus, 
and algebraic topology as its discrete counterpart. 

The use of differential geometry and algebraic topology enjoys increasing popularity for the 
development of compatible schemes, [5, 6, 11, 12, 13, 21, 34, 35]. Compatible discretizations are 
often combined with mixed formulations. Mixed formulations are described extensively in among 
others [15, 30] and in terms of differential forms in [5, 6] for the Hodge-Laplacian and in [40] for 
the VVP formulation of the Stokes problem. 

The MMSEM contains compatible finite elements that are compatible with all admissible types 
of boundary conditions for the Stokes problem in VVP formulation. We will show that the 
method obtains optimal rates of convergence for all variables on curvilinear meshes and for all 
admissible boundary conditions, i.e. standard and nonstandard. It is therefore extending the error 
estimates found in literature, which are often specifically constructed for certain types of boundary 
conditions, [1, 4, 10, 14, 24, 29]. To show optimal convergence a priori error estimates are derived. 

This is an improvement with respect to the well-known Raviart-Thomas compatible finite ele- 
ments. These are not compatible in case of Dirichlet boundary conditions and therefore lead to 
suboptimal convergence behavior, as was shown in [4, 24]. This non-compatibility results in a 
decrease in rate of convergence of maximal | order. 

From a physical/fluid dynamics point-of-view the new method is relevant because it combines 
optimal convergence with a pointwise divergence- free discretization (in absence of any mass source) 
of arbitrary order on curvilinear meshes, valid for all allowable types of boundary conditions, among 
which the no-slip condition. 

The derived rates of convergence are confirmed using simple manufactured solution problems, 
discretized on both Cartesian and curvilinear meshes. The fact that the analysis holds for all 
admissible boundary conditions is also reflected in the numerical results. 

This paper is organized as follows: First an introduction into differential geometry is given 
and the Stokes problem is reformulated in terms of differential forms. In Section 3 the mixed 
formulation is given and well-posedness is proven. In Section 4 the key properties of the mimetic 
discretization arc explained that lead to compatible function spaces. This includes a discussion on 
the relevant properties of algebraic topology, the definitions op mimetic operators, the introduction 
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of mimetic spectral element basisfunctions and finally the proof of discrete well-posedness. Having 
formulated the conforming/compatible finite element spaces, the error estimates are developed in 
Section 5 and the numerical results are shown in Section 6. 

2. Notation and preliminaries 

2.1. Differential forms. Differential forms offer significant benefits in the construction of structure- 
preserving spatial discretizations. For example, the coordinate-free action of the exterior derivative 
and generalized Stokes theorem give rise to commuting properties with respect to mappings be- 
tween different manifolds. Acknowledging and respecting these kind of commuting properties are 
essential for the structure preserving behavior of the mimetic method. 

Only those concepts from differential geometry which play a role in the remainder of this paper 
will be explained. More can be found in [2, 26, 27, 41]. 

Let A k (Sl) denote a space of differential k-forms or k-forms, on a sufficiently smooth bounded 
n-dimensional oriented manifold f2 C K™ with boundary T = dSl. Every element a e A k (fl) has a 
unique representation of the form 

(2.1) a = //(x)da; Jl A dx 12 A • • • A dx ik , 

where I = i\, . . . with 1 < i\ < . . . < if. < n and where //(x) is a continuously diffcrentiable 
scalar function, //(x) e C°°(f2). Differential fc-forms are naturally integrated over fc-dimensional 
manifolds, i.e. for a € A k (il) and flk C K™, with k = dim(Slfc), 

(2.2) / ael (a,fi fc )eR, 

where (•, •) indicates a duality pairing between the differential form and the geometry. Note that 
the n-dimensional computational domain is indicated as fl, so without subscript. The differential 
forms live on manifolds and transform under the action of mappings. Let $ : — » f2 be a 
mapping between two manifolds. Then we can define the pullback operator, $* : A k (n) — > A k (il), 
expressing the fc-form on the n-dimcnsional reference manifold, 0. The mapping, and the 
pullback, $*, are each others formal adjoints with respect to a duality pairing (2.2), 

(2.3) [ a= [ $*a & (a,$(fij)> = <**<*» "l>> 

where is an Z-dimensional submanifold of 17 and = <I> (O^ ) a fc-dimensional submanifold of 
£1 A special case of the pullback operator is the trace operator. The trace of fc-forms to the 
boundary, tr : A k (fl) — > A fe (<9£l), is the pullback of the inclusion of the boundary of a manifold, 
dn 0, see [41]. 

The wedge product, A, of two differential forms a e A fe (fi) and b £ A l (Q) is a mapping: 
A : A fe (f2) x A'(fi) -> A fe+i (£l), k + I < n. The wedge product is a skew-symmetric operator, i.e. 
aAb= (-l) kl bAa. 

An important operator in differential geometry is the exterior derivative, d : A fc (£l) — > A fe+1 (f2). 
It is induced by the generalized Stokes' theorem, combining the classical Newton-Leibnitz, Stokes 
circulation and Gauss divergence theorems. Let ftk+i be a (fc + l)-dimensional manifold and 
a e A fe (f2), then 

(2.4) / a= da ^ (a, dQ, k+1 ) = (da, fi fe+ i), 

J dflk + i JCik+i 

where dflk+i is a fc-dimensional manifold being the boundary of fifc+i. The duality pairing in (2.4) 
shows that the exterior derivative is the formal adjoint of the boundary operator d : Slk+i ~^ ^fe- 
The exterior derivative is independent of any metric and coordinate system. Applying the exterior 
derivative twice always leads to the null (fc + 2)-form, d(da) = 0( fe+2 ' for all a e A fc (fi). As a 
consequence, on contractible domains the exterior derivative gives rise to an exact sequence, called 
De Rham complex [27], and indicated by (A, d), 

(2.5) R ^ A°(fi) -A A : (0) 4A"(!])Ao. 
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In vector calculus a similar sequence exists, where, from left to right for R 3 , the d's denote 
the vector operators grad, curl and div. The exterior derivative and wedge product are related 
according to Leibnitz's rule as: for all a <G A k (il) and b <G A'(f2), 

(2.6) d(oAi)) = daA&+(-l)' £ aAd6, for k + I < n. 

The pullback operator and exterior derivative possess the following commuting property, 

(2.7) $ i da = d$ 4 o, VaeA fe (0). 

In this paper we will consider Hilbert spaces L 2 A k (Q) D A fe (f2), where in (2.1) the functions 
//(x) e L 2 (Q). The pointwise inner-product of fc-forms, (•, •), is constructed using inner products 
of one- forms, that is based on the inner product on vector spaces, see [26, 27]. The wedge product 
and inner product induce the Hodge-* operator, * : L 2 A k (VL) — > L 2 A™~ fe (0), a metric operator 
that includes orientation. Let a, b e L 2 A k (fl), then 

(2.8) a A*b := (a,b)a, 

where a e A"(Q) is a unit volume form, u = *1. In geometric physics the Hodge-* switches 
between an inner-oriented description of physical variables and an outer-oriented description. 
See [40, 41, 43, 52] for a thorough discussion on the concepts of inner and outer orientation. 
The space of square integrable fc-forms on can be equipped with a L 2 inner product, (•, -) n : 
L 2 A fe (0) x L 2 A k (n) -> R, given by, 

(2.9) K & )n := / (a,b)a^ = [ aA*o. 

Jn Jn 

The norm corresponding to the space L 2 A k (Q) is ||a||x,2A fc = \J ( a i a ) n - Higher degree Sobolev 

spaces, H m A k , consists of all fc-forms as in (2.1) where //(x) e H m (Q), with corresponding norms 
\a\ HmA k and || The Hilbert space associated to the exterior derivative HA k (il) is defined 

as 

(2.10) HA k (n) = {ae L 2 A k {n) | da e L 2 A fe+1 (0)}. 

and the norm corresponding to HA k (il) is defined as ||a|| 2 ^ Afc := ||a|| 2 2 ^ k + ||da|| 2 2A fc + i • The HA k - 
semi-norm is the L 2 -norm of the exterior derivative, |a|#A fc = ||da||z,2 A fc+i . Note that H 1 A k (Q) C 
HA k (fl) C L 2 A k (fl), where the left equality holds for k = and the right for k = n. The L 2 -de 
Rham complex, also called Hilbert complex [16], (HA, d), is the exact sequence of maps and spaces 
given by 

(2.11) R ^ HA°(Q) A HA\n) A ±+ HA n (Q) A 0. 

In terms of vector operations the Hilbert complex becomes for ft C R 3 , 

H\n) ^ if (curl, ft) ^4 (div.fi) ^> £ 2 (ft), 

and for del 2 , either 

H\n) ^ (rot, ft) ^ L 2 (ft), or if^ft) ^4 (curl, ft) L 2 (ft). 
The two are related by the Hodge-* operator (2.8), see [46], 

FA°(ft) A if A 1 (ft) A i 2 A 2 (ft) if^ft) ^4 if(curl.ft) L 2 (ft) 

(2.12) *^ *^ *^ *i *j; ★ $ 

L 2 A 2 (ft) A ffA^ft) A ifA°(ft) i 2 (ft) ^ if (rot, ft) 1^ if x (ft). 

Remark 1. The upper complex is associated with outer- oriented k-forms, i.e. k-forms that are 
associated with outer- oriented manifolds, and the lower complex is associated with inner-oriented 
k-forms. In this paper we mainly consider the upper complex and circumvent the lower complex 
by means of integration by parts. Only the pressure and tangential velocity boundary conditions 
are given on the lower complex, as we will see in the following sections. 
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A similar double Hilbert complex can be constructed in R 3 . Since the exterior derivative is 
nilpotcnt, it ensures that the range, B k := di7A fc_1 (ri), of the exterior derivative on (k— l)-forms 
is contained in the nullspace, Z k := { a e HA k (Q) \ da = }, of the exterior derivative on 
fc-forms, B k CZ k . 

Every space of fc-forms in the complex (HA, d) can be decomposed into the nullspace of d, and 
its orthogonal complement, HA k (fl) = Z k © Z k ' ± . This is the Hodge decomposition, where on 
contractible domains Z k = B k . By the Hodge decomposition it follows that the exterior derivative 
is an isomorphism d : Z k '^ — »■ B k+1 . 

The inner product gives rise to the formal Hilbert adjoint of the exterior derivative, the codiffcr- 
ential operator, d* : H*A k (n) -> L 2 A fc_1 (fl). Let H*A k (n) = {«£ L 2 A k (n) \d*a e L 2 A k ~ 1 (Ct)}, 
then 

(2.13) (da,b) n = (a,d*b) Q , Va e #A fc-1 (fi), b e H*A k . 

In case of non-zero trace, and by combining (2.9), (2.4) and (2.6), we obtain integration by parts, 

(2.14) (a,d*b) n = (da,b) n - traAtr*6. 

Jan 

Also the codiffcrcntial operator is nilpotent, d*(d*a^^) = 0, i.e., its range is contained in its 
nullspace, B*- k C Z*- k , where B*' k := d*H*A k+1 (Q) and Z*' k := { a e H*A k (n) \ d*a = }. 
In fact the codifferential is an isomorphism d* : Z*' fe,J - — > B*' k+1 , where Z*- k,± - follows from the 
following Hodge decomposition, A k (Q) = Z*' k © Z*^' 1 - . On contractible manifolds this gives rise 
to the following exact sequence, 

(2.15) H*A°{Q) i?*A x (0) <^ A H*A n ({l) ^ R. 

In vector notation from right to left the d*'s denote the -grad*, curl* and -div* operators in 
R 3 , as were also mentioned in the introduction. However, whereas the exterior derivative is a 
metric-free operator, the codifferential operator is metric-dependent. The Hodge-Laplace operator, 
A : H 2 A k (fi) — > L 2 A k (VL), is constructed as a composition of the exterior derivative and the 
codifferential operator, 

(2.16) - A a := (d*d + dd*)a. 

An important inequality in stability analysis, relating the L 2 A fc -norm and the HA k -norm, is 
Poincare inequality. 

Lemma 1 (Poincare inequality). [6] Consider the Hilbert complex (HA,d), then the exterior 
derivative is a bounded bijection from Z k - ± to B k+1 , and hence, by Banach's bounded inverse 
theorem, there exists a constant cp such that 

(2.17) ||a|| HAfc <cp||do||i2 A * + i, VaeZ^. 

Finally, for Hilbert spaces with essential boundary conditions we write, H A k (ft) := {a € 
HA k (Cl) I tr a — 0}, and for natural boundary conditions we consider the following trace map, 
tr * : HA k (Q) -> H^A n - k {dQ). 

2.2. Stokes problem in differential form notation. Consider again a bounded contractible 
domain fid". Because we require exact conservation of mass and because we can perform exact 
discretization of the exterior derivative, see Section 4.2, we use the following formulation for the 
Stokes problem: let (w, u,p) G {A™- 2 (0) x A"" 1 ^) x A™(0)}, then the VVP formulation is given 



by 

(2.18a) oj-d*u = Q, inA"- 2 (ft), 

(2.18b) d*du + dw + d*p = f, inA™" 1 ^), 

(2.18c) du = g, inA n (ft). 



In the VVP formulation the pressure in (2.18b) acts as a Lagrange multiplier for the constraint 
on velocity, (2.18c), whereas velocity in (2.18a) acts as a Lagrange multiplier for the constraint on 
vorticity in (2.18b). 
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Let T = dSl be the boundary of il, where 

r = r u ur ( , r w nr t = 0, and r = r„ur 7r , r„nr, = 0. 

We will impose the tangential vorticity and normal velocity as essential boundary conditions, 
and the tangential velocity and the pressure plus divergence of velocity as the natural boundary 
conditions: 



(2.19a) 


tr u> 


= 


on r w , 




(2.19b) 


tr u 


= 


on r„, 




(2.19c) 


tr *u 


= u b , t 


on r t , 


with u b . t e H^A\T t ), 


(2.19d) 


tr * (du + p) 


= n 6 


on 1^, 


with n fc e H^A°(r n ). 



Then the boundary T can be partitioned into four sections, T = {J i=1 Ti, with n Tj =0 for 
i 7^ j, where 

(2.20) r i: =r t = r„, r 2 :=r t = r T , r 3 :=r w = r„, r 4 :=r w = rv 

This decomposition, introduced before in [23, 36, 40], shows all admissible boundary conditions. 
It will also follow directly from the mixed formulation, see (3.7), Section 3. 

In case, T = Ti n T 3 , T 2 U T 4 = 0, no pressure boundary conditions are prescribed, and so the 
pressure is only determined up to an element p e Z*' n , i.e. up to a constant. As a post processing 
step either the pressure in a point in SI can be set, or a zero average pressure can be imposed; i.e. 
J Q p = 0. In case T = T 4 , no velocity boundary conditions are prescribed, and so the solution of 
velocity is determined modulo a curl*-free element, i.e. modulo u e 

3. Mixed formulation 

3.1. Mixed formulation of Stokes problem. The use of a mixed formulation is based on the 
following reasoning; We know how to discretize exactly the metric-free exterior derivative d, but 
it is less obvious how to treat the codiffcrcntial operator d* . 

3.1.1. Generalized Poisson problem. Take for example the generalized Poisson problem using the 
Hodge-Laplacian acting on fc-forms, (dd* + d*d)w = /, on Q with boundary T = dtt. A standard 
Galerkin approach, using integration by parts (2.14), would give; find u e HA k (Q) n H*A k (n) 
with du G H*A k+1 (Q) and d*u e HA k ~ 1 (n), given / e L 2 A k (Vt), such that 

(3.1) (d*v,d*u) n + (dv,du) Q = (v,f) n , Vv e HA k (n) n H*A k (n). 

It has a corresponding minimization problem for an energy functional over the space HA k (Sl) (~l 
H*A k (fl). The standard Galerkin formulation is coercive, which immediately implies stability. 
Corresponding to this standard Galerkin formulation one usually chooses a 7? 1 A fe (rj)-conforming 
approximation space. This could be a standard continuous piecewise polynomial vector space 
based on nodal interpolation. 

However, in case of a nonconvex polyhedral or curvilinear or noncontractiblc domain fl, for 
allmost all /, HA k (Q,) n H*A k {n) <£ H 1 A k (Q,). Consequently, the solution will be stable but 
inconsistent in general, [20]. In other words, the solution converges to the wrong solution. Unfor- 
tunately, it seems not possible to construct HA k (Sl) n H*A k (Ct) conforming finite element spaces. 
Alternatively, one proposed to use mixed formulations, [15]. In contrast to standard Galerkin, the 
mixed formulation uses integration by parts (2.14) to express each codifferential in terms of an 
exterior derivative and suitable boundary conditions. 

Consequently, mixed formulations require only HA k (O)-conforming finite element spaces, which 
are much easier to construct. Therefore, in all cases mixed formulations do converge to the 
true solution. Mixed formulations correspond to saddle point problems instead of minimization 
problems. 

The derivation of the mixed formulation of the Poisson problem consists of three steps: 
(1) Introduce an auxiliary variable w = d*u in _ffA fe_1 , 
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(2) multiply both equations by test functions (r, v) <E {HA k ~ 1 x HA k } using L 2 -inner prod- 
ucts, 

(3) use integration by parts, as in (2.14), to express the remaining codiffcrcntials in terms of 
the exterior derivatives and boundary integrals. 

Again the boundary may constitute up to four different types of boundary conditions, 

(3.2a) trw = on T u , 

(3.2b) tr u = on T„, 

(3.2c) tT*u = u b>t on r f , with u M € H^A n ~ k (T t ), 

(3.2d) tr*du = g b on IV, with g b e i?' A"~ fe_1 (iV). 

Then also for the generalized Poisson problem the boundary T consists up to four sections as 
defined in (2.20). To obtain a unique solution for the corresponding mized formulation, we define 
the following Hilbert spaces, 

(3.3) W:={TeiTA fc - 1 (fi)|trT = Oonr ai }, 

' {v e HA k (n) | tr v = on T„ } if T 4 = 0, 
{ v G HA k (n)\Z*< k (n) | tr v = on r„ } if r 4 7^ 0, 



(3.4) V := • 



with corresponding norms, || • \\w, \\ ■ \\v, respectively. The resulting mixed formulation for the 
Poisson problem for all < k < n becomes: find (cj,u) € {W x V}, given / £ L 2 A k , for all 
(t,v) e {W x V}, such that 

(3.5a) (r,w) f2 - (dr,u) f2 = - / trrA« w , 

•Mura 

(3.5b) (v,dw) Q + (dv,du) n = (v,f) n + [ tivAg b . 

Jr 2 ur 4 

Note that, for a scalar Poisson, it is not a choice whether to use Galerkin or mixed formulation, 
but it depends on whether the scalar is a 0-form or an n-form. This is determined by the physics. 

3.1.2. Stokes problem. In a similar way the mixed formulation of the VVP formulation of the 
Stokes problem is obtained. Consider the Hilbert spaces W and V defined in the previous section, 
where k — n — 1, and define the following Hilbert space 



(3.6) Q :-- 



q eL 2 A n {fi), ifr w ^0, 
q e L 2 A n {n)\z*' n , if r T = 0, 



with corresponding norm j| • \\q and where Z*' n = R. Then the mixed formulation of the VVP 
formulation reads: find (uj, u,p) e {W x V x Q}, for the given data / e L 2 A Il_1 (fi), g e L 2 A"(fi) 
and natural boundary conditions u btt € A n_1 (ri), n fc e H2A n (T n ), for all (r, v,q) e {W x 
V x Q}, such that 

(3.7a) (r,u) n - (dr,u) n = - trrAu M , 

Jriur 2 

(3.7b) (v,duj) n + (dv,du) n + (dv,p) n = («,/)„+ / trwAn fc , 

Jr 2 ur 4 

(3.7c) (g,d«) n =(g,</) n . 
Again use is made of integration by parts, (2.14). 

Proposition 1. [7] Problems (2.18)-(2.19) and (3.7) are equivalent, in the sense that any triple 
(to,u,p) € {W x V x Q} is a solution of problem (2.18)-(2.19) if and only if it is a solution of 
problem (3.7). 
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3.2. Well-posedness of mixed formulation. Before we continue first define the following 
nullspaces of W, 

(3.8a) Z w := { r e W | dr = }, 

(3.8b) Z* w := { t e W | d*r = }, 

and consider the following decompositions, W = Z^y © Z^ and W = Z^ © Zt^~ . Since vorticity 
is defined as uj = d*u, we have uj € ZJy, and because we consider contractiblc domains only, it 
follows that uj e Z^r. Note that for n = 2, Z{y = W. A similar decomposition can be made for 
V. Define 

(3.8c) Z v := { v e V | dv = }, 

then V = Zy © Zy . The velocity is decomposed as u = uz + u± , where 112 € Zy and u± e Zy . 

We can write the mixed formulation of (3.7) in a more general representation, using four con- 
tinuous bilinear forms, 

a(-, ■):=(■, On : W x W -»■ R, b(-, •) := (d, -)n ^xQ^R, 
c(-,-):=(d-,-)n :^x|/^K, e(-, •) := (d, d-)n : V x V -+ M, 
and three continuous linear forms 

f(-):=0./)n+ / tr -An b :F^M, 
Jr 2 ur 4 

g(-) := (■, 9)o : Q ->• R, h(') := - / tr • Au M : W ->■ K. 

Jriur 2 

The mixed formulation becomes 

(3.9a) a(r, uj) - c(r, u) = h(r), Vr e W, 

(3.9b) e(u,u) +c(w,u) + b(u,p) = f(v), Vv e V, 

(3.9c) b(u,«)=g(«). VqeQ. 

There exists continuity constants < c a , Cb, c c , c e < 00 such that 
(3.10) 

a(r, k) < c b \\t\\w\\k\\w, g) < Cb||v||v||g||Q, c(t, u) < c c ||r||w||w|jy, e(v, w) < c e \\v\\y\\w\\y. 

By Cauchy-Schwarz we know that c a = 1, however we write c a for generality purpose. The 
continuous linear forms are bounded such that 

(3-11) f(v)<\\f\\\\v\\ v , g(v)<\\g\\\\v\\ v , h(^) < II^HII-rllwr. 

At first restrict to all v = v z e Z v . This gives the vorticity- velocity subproblem, which is a saddle 
point problem: 

(3.12a) a(r,w) -c(t,u z ) = h(r), Vr e W, 

(3.12b) c(w z ,cj) = f(u^), Vw z e Zy. 

Proposition 2. [23] System (3.12) has a unique solution (uj,Uz) € {W^ x Zy} if there exists 
positive constants a, 7, suc/i that we have coercivity in the kernel of W, 

(3.13) inf sup - n — n— > a, inf sup - n — n — > a, 

T Z eZ w KzeZw \\Tz\\w\\KZ\\W kz£Z w TZ £Z W \\Tz\\w\\KZ\\W 

and satisfies the following inf-sup condition for c(t,v z ), 

(3.14) inf sup C } T \ VZ \, > 7, 

v z ez v T £-\y \\t\\w\\ v z\\v 
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Proof. The proof of (3.13) is straightforward, see e.g. [15]. For (3.14), we have c(t,v z ) = 
(dr, Vz)n, where d : — > Zy. Thus, given vz £ Zy there exists a unique t v £ Z^ such 
that dr v = vz and ||t„||v(/ < cp||u^||y by Lemma 1. Therefore 

c(r,v z ) c(t v ,v z ) \\vz\\v . 1 n n 

SUp ——r > — jj = j j— - > H^llv- 

tEW \\ T \\W \\ T v\\w \\ T v\\w C P 

□ 

Proposition 3. The full problem (3.9) has a unique solution (u,u,p) £ {W, V, Q} if conditions 
(3.13) and (3.14) from Proposition 2 are satisfied and additionally is there exists positive constants 
e, (3, such that we have coercivity in the range of V, 

• t e(vj_,wj_) . e(v±,w±) 

(3.15) ml sup n — n n — > e, ml sup j r r r — > £, 

v±ez^ w±ez ± \\vj_\\ v \\wj_\\ v w± ez$ v±eZ .L \\vj_\\ v \\wj_\\ v 

and satisfies the following inf-sup condition for b(v, q), 

(3.16) inf sup „ b j V ' q \ >(3>0. 

qeQ veV \\v\\v\\q\\Q 

Proof. The proof is similar to that of Proposition 2. See also [8], Section 7.1. □ 

So well-posedness of the Stokes problem (3.7) relies only on the Hodge decomposition and the 
Poincare inequality. 

Corollary 1. [7, 23] Problem (3.7) is well-posed according to Propositions 2 and 3. That is, for any 
given data f £ L 2 A" _1 (fi) and g £ L 2 A n (fl) and natural boundary conditions ut,j € A n ~ 1 (T t ) 
and lib £ if 5 A" (IV) , there exists a unique solution (oj,u,p) £ W x V x Q satisfying (3.7). 
Moreover, this solution satisfies: 

(3.17) \\u\\ w + \\u\\ v + HpHo < C (||/||l»a»-i + HfflU»A» + ||«Mll ff J + W\ H \ Kn ) , 
where C is a constant depending only on the Poincare constant cp and the continuity constants. 

4. Compatible spectral discretization 

Well-posedness of the Stokes problem in VVP formulation relies solely on the Hodge decom- 
position and the Poincare inequality. For a compatible discretization, these properties need to be 
respected as well in the finite dimensional spaces. Key ingredient to obtain a discrete Hodge de- 
composition and discrete Poincare inequality is the construction of a bounded projection operator 
that commutes with the exterior derivative. 

The compatible spectral discretization consists of three parts. First, the discrete structure 
is described in terms of chains and cochains from algebraic topology, the discrete counterpart 
of differential geometry. This discrete structure mimics many of the properties from differential 
geometry. Secondly, mimetic operators are introduced that relate the continuous formulation in 
terms of differential forms to the discrete representation based on cochains and finite dimensional 
differential forms. Thirdly, mimetic spectral element basis functions are described following the 
definitions of the mimetic operators. In this paper we address these topics only briefly. More 
details of the mimetic spectral element method can be found in [40, 41]. Finally, well-posedness 
of the discrete numerical formulation is proven and interpolation error estimates are given. 

4.1. Algebraic Topology. Let D be an oriented cell-complex covering the manifold O, describing 
the topology of the mesh, and consisting of fc-cells t^), fc = 0, . . . ,n. The two most popular classes 
of fc-cells in literature to describe the topology of a manifold are either in terms of simplices, see 
for instance [44, 51, 53], or in terms of cubes, see [42, 52]. From a topological point of view both 
descriptions are equivalent, see [22]. Despite this equivalence of simplicial complexes and cubical 
complexes, the reconstruction maps in terms of basis functions, to be discussed in Section 4.2, 
differ significantly. For mimetic methods based on simplices see [5, 6, 21, 47], whereas for mimetic 
methods based on singular cubes see [3, 37, 39, 50]. We restrict ourselves to /c-cubes, although we 
will keep calling them fc-cells. 
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The ordered collection of all fc-cells in D generate a basis for the space of fc-chains, Cfc(-D). 
Then a fc-chain, c^) e Cfc(-D), is a formal linear combination of fc-cells, t^),% € D, 

(4.1) c( fe ) = ~y]ciT(k),i- 

i 

The boundary operator on fc-chains, d : Ck(D) — > Ck-i{D), is an homomorphism defined by 
[33, 44], 

(4.2) dc(k) = dy]ciT(fc),i : = ^ Cl <9 (r^) . 

i i 

The boundary of a fc-cell T( fc ) will then be a (fc — l)-chain formed by the oriented faces of r^y 
Like the exterior derivative, applying the boundary operator twice on a fc-chain gives the null 
(fc — 2)-chain, ddc^ = 0( fe _ 2 ) for all C( fc ) e Ck{D). The set of fc-chains and boundary operators 
gives rise to an exact sequence, the chain complex (Ck(D),d), 

( 4 -3) ••• ^— C fc _i(D) Ck(D) C fc+ i(D) ••• . 

Let i3fe be the range and Z k be the nullspace of d in Cfe. Then the topological Hodge decomposition 
of the space of fc-chains is given by = ® Zjjj-, where Z k = B k on contractible domains 1 . The 
boundary operator on chains in (4.3) is a bijection that maps d : Z k — > Bk-i- 

Dual to the space of fc-chains, Ck(D), is the space of k-cochains, C k (D), defined as the set 
of all linear functionals, c^' : Ck(D) — > R. The duality is expressed using the duality pairing 
(c( k \ C(fc)) := c( fe ^(c( fe )). Note the resemblance between this duality pairing and the integration of 
differential forms (2.2). 

Let {T(k),i\ form a basis of Ck(D), then there is a dual basis {r^' 2 } of C k (D), such that 
T ^' l ( T (k),i) = <5) and all fc-cochains can be represented as linear combinations of the basis elements, 

(4.4) c<*> = ^cM fe K 

With the duality relation between chains and cochains, we can define the formal adjoint of the 
boundary operator which constitutes an exact sequence on the spaces of fc-cochains in the cell 
complex. This formal adjoint is called the coboundary operator, 8 : C k (D) — > C k+1 (D), and is 
defined analogous to (2.4) as 

(4.5) (Sc^,c (k+1) ) := (c«, dc (fc+1) ), Vc« e C k (D) and Vc (fe+1) e C k+1 (D) . 

Note that expression (4.5) is nothing but a discrete Stokes' theorem and that the coboundary 
operator is nothing but a discrete exterior derivative. Also the coboundary operator satisfies 
88c^ = 0( k+2 \ for all c( fc ) <E C k (D), and gives rise to an exact sequence, called the cochain 
complex (C k (D),8), 

(4.6) ... C k -\D) — i-> C k {D) — i-> C k+1 (D) — 

Let B k be the range and Z k be the nullspace of 8 in C fe , then a Hodge decomposition of the 
space of fc-cochains is given by C k — Z k Z k ' ± , where Z k = B k on contractible domains. The 
coboundary operator in (4.6) is a bijection that maps 8 : Z k '^ — > B k+1 . Note the similarity 
between this map, that of the boundary operator on fc-chains and that of the exterior derivative 
on fc- forms. 

4.2. Mimetic Operators. The discretization of the flow variables involves a bounded projection 
operator, Wh, from the complete space HA k (Q) to a conforming subspace A^(f2;Cfc) C HA k (£l). 
The projection operation consists of two steps, a reduction operator, 1Z : HA k (fl) — > C k (D), 
that integrates the fc-forms on fc-chains to get fc-cochains, and a reconstruction operator, X : 
C k (D) — > A k (Q;Ck), to reconstruct fc-forms from fc-cochains using appropriate basis- functions. 



Although 'perpendicular' in a topological space is not well defined, we refer to Zjjj- as the complement space of 
Zk m Cfe- 
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These mimetic operators were already introduced before in [11, 38]. A composition of the two 
gives the projection operator 7r^ = X o TZ as is illustrated below. 

HA k (n) A k (Q;C k ) 
TZ 




C k (D) 

These three operators together constitute the mimetic framework. An extensive discussion on 
mimetic operators can be found in [40, 41]. 

The reduction TZ and reconstruction X operators are defined below. The fundamental property 
of TZ and X is the commutation with differentiation in terms of exterior derivative and coboundary 
operator. 

The reduction operator TZ : HA k (Q) — > C k (D) is a homomorphism that maps differential forms 
to cochains. This map is defined by integration as 

(4.7) (TZa,T (k) ) := [ a, Va G HA k (Q,), T( fe ) G C k (D). 

Then for all C( fe ) G C k (D), the reduction of the fc-form, a G HA k (fl), to the fc-cochain, a( fe ) G 
C k (D), is given by 

(4.8) a< fc >(c (fc) ) := {Ka,c (k) ) ( = } £ c^TZa, r (fe);i ) ( = 7) ^ f a = f a. 

The reduction maps has a commuting property with respect to differentiation in terms of 
exterior derivative and coboundary operator, 

(4.9) TZd = 5TZ, onHA k {Vi). 

Since TZ is defined by integration, (4.9) follows directly from Stokes theorem (2.4) and the duality 
property (4.5). 

Next by definition also the reconstruction map X : C k (D) — > A^(f2;Cfc) needs to have a com- 
muting property with respect to differentiation in terms of exterior derivative and coboundary 
operator, 

(4.10) dX = XS, on C k (D). 

The reconstruction X must be the right inverse of TZ, so TZX — Id on C k (D), and we want it to be 
an approximate left inverse of TZ, so XTZ = Id + 0(h p ) on HA k (fl). This composition is defined 
as the projection operator. 

Definition 1 (Bounded projection operator). The composition XoTZ will denote the projection 
operator, ir^ := XTZ : HA k (Q) — > A^(0; C k ), allowing for a finite dimensional representation of a 
k-form, 

(4.11) 7r h a := XTZa, ir h a G A k h {tl; C k ) C HA k (n). 

where XTZa is expressed as a combination of k-cochains and interpolating k-forms. The projection 
operator ir h is a bounded operator if for C < oo and for all a G A k (il) we have ||7r^a|| ffA fe < 
C\\a\\ HA k. 

A proof that Wh is indeed a projection operator is given in [41]. In Section 4.5 also boundedness 
is proven. 

Lemma 2 (Commutation property). There exists a commuting property for the projection 
and the exterior derivative, such that 

(4.12) dTr h = ir h d onHA k ({l). 
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Proof. Express the projection in terms of the reduction and reconstruction operator, then 
dn h a (4 = ] dZTLa ( " 0) 1611a ( = 9) Tilda (4 = 1} ir h da, Va e HA k (n). 

□ 

Note that it is the intermediate step 1511a that is used in practice for the discretization, sec 
[40], and (4.20) on page 14. 

Corollary 2 (Discrete Hodge decomposition). From Lemma 2 it follows that B^ :— nf l B k C 
B k , Zfi := TT h Z k C Z k and that on contractible domains, Z^ := ■KhZ k ' ± C Z k < ± . Then the 
discrete Hodge decomposition is given by A k — Z k (B Z^'' 1 . As a consequence of Lemma 2 and the 
discrete Hodge decomposition we have Z\y h C Z\y , Zy h C Zy , dWh C Vh and dVh = Qh, which 
shows that the discretization method is compatible. 

Finally, we do not restrict ourselves to affine mappings only, as is required in many other 
compatible finite elements, like Nedelec and Raviart-Thomas elements and their generalizations 
[5, 45, 48], but also allow non- affine maps such as curvilinear transfinite or isoparametric mappings 
of quadrilaterals or hexahcdrals, [31], where $ and its inverse are piecewise sufficiently smooth, 
i.e. 

(1) $ is a C p+1 -diffcomorphism, 

(2) \4>\ WL <Ch l , l<p+l, 

(3) I^Vi, < Ch~\ l<p + l. 

This allows for better approximations in complex domains with curved boundaries, without the 
need for excessive refinement, while maintaining design convergence rates, [19]. This is possible 
since the projection operator irh commutes with the pullback <f>*, 

(4.13) <f>*ir h = ?r h $* on HA k (Q). 

An extensive proof is given in [41]. 

4.3. Numerical stability. Essential ingredients in proving numerical stability are the discrete 
Hodge decomposition and the discrete Poincare inequality. Because the complexes (HA,d) and 
(A^,d) are each others supercomplex and subcomplex, respectively, the discrete Poincare inequal- 
ity is directly related to the Poincare inequality in Lemma 1 and the bounded projection in 
Definition 1. 

Lemma 3 (Discrete Poincare inequality). Let (HA,d) be a bounded closed Hilbert complex, 
(A^,d) a subcomplex, and ir^ a bounded projection. Then 

^ 4 I IKHifA* < c Ph \\da h \\ L 2 A k, a h <^Z k ^, 

\ 1 < Cph < c P . 

Proof. Given ah € Z h ' . From the Hodge decomposition and the bounded projection it follows 
that B k h C B k and Z£ C Z k and from the commutation relation (4.12) it follows that Z^ C Z k ^ . 
Since we consider a proper subspace, Lemma 1 is still valid, with cpy l < cp. □ 

Theorem 1 (Discrete well-posedness). Let (A^,d) be a subcomplex of the closed Hilbert 
complex (HA,d). Then there exists constants cth, Ph^h> depending only on Cph, such that for 
any (rh,Vh,qh) € Wh x Vh x Qh, there exists a stable finite dimensional solution (ujh,Uh,Ph) € 
WhX Vh x Q h of the Stokes problem (3.7), with 

(4.15) a h > a > 0, (3 h > P > 0, 7^ > 7 > 0. 

Proof. This is just Propositions 2 and 3 applied to the complex (A^,d), combined with the fact 
that the constant in the Poincare inequality for A^ is cph < cp by Lemma 3. □ 
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4.4. Mimetic spectral element basis-functions. The finite dimensional differential forms used 
in this paper are polynomials, based on the idea of spectral element methods, [18]. The mimetic 
spectral elements used here were derived independently in [28, 49], and are more extensively 
discussed in [41] . Only the most important properties of the mimetic spectral element method are 
presented here. 

In spectral element methods the domain Q is decomposed into M non-overlapping, in this case 
curvilinear quadrilateral or hexahedral, closed sub-domains Q m , 

M 



^ = U Qm, Qm n Ql = 8Q m fl 8Ql, m/1, 



m— 1 

where in each sub-domain a Gauss-Lobatto grid is constructed. The complete mesh is indicated 

by 2 : = Em=l<5m- 

The collection of Gauss-Lobatto meshes in all elements Q m € Q constitutes the cell complex 
D. For each element Q m there exists a sub cell complex, D m . Note that D m fl Di, m ^ I, is not 
an empty set in case they are neighboring elements, but contains all fc-cells, k < n, of the common 
boundary. 

Each sub-domain Q m € Q is mapped from the reference element, Q = [—1, 1]™, using the 
mapping $ m : Q — > Q m . Then all flow variables defined on Q m are pulled back onto this reference 
element using the following pullback operation, <j>^ : A^(Q m ;Ck) — > Ajj(Q;Cfe). 

The basis-functions that interpolate the cochains on the quadrilateral or hexahedral elements 
are constructed using tensor products. It is therefore sufficient to derive interpolation functions 
in one dimension and use tensor products afterwards to construct n-dimensional basis functions. 
A similar approach was taken in [17]. Because projection operator and pullback operator com- 
mute (4.13), the interpolation functions are discussed for the reference element only. Since the 
mappings Q m and their inverse are assumed to be sufficiently smooth, the rates of convergence for 
interpolation estimates on the physical elements are equal to that of the reference element. Only 
the constants C that will appear below will depend on the mappings <J> m , but will be independent 
of the meshsize and polynomial order. 

Consider a 0-form a e HA°(Q) on Q := £ e [—1, 1], on which a cell complex D is defined that 
consists of N + 1 nodes, where — 1 < £ < . . . < <^v < 1, an d N edges, r^)^ = [&_!,&], of 
which the nodes are their boundaries. Corresponding to this set of nodes (0-chains) there exists a 
projection using N th order Lagrange polynomials, k(£,), to approximate a 0-form, as 

N 



(4.16) 7T h a = J2 a ih(0- 



j=0 



The property of Lagrange polynomials is that they interpolate nodal values. They are therefore 
suitable to reconstruct a 0-form form the 0-cochain a^ ' = IZa, a e A°(f2), containing the set 
a i = a (&) for i = 0, . . . , iV. Lagrange polynomials are in fact 0-forms, ?»(£) € A° (Q; Co). Lagrange 
polynomials are constructed such that their value is one in the corresponding point and zero in 
all other grid points, 



(4.17) TZh(0 = Ufa) 



1 if i = p 
if i ^ p 



In [28, 49] a similar basis for projection of 1-forms was derived, consisting of 1-cochains and 1- 
form polynomials, that is called the edge polynomial, ej(£) G A^(Q). Let b e L 2 A 1 (Q), then the 
projected 1-form is given by 



N 



(4-18) 7r fc &(0 = I> e *(0. 



1=1 
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where the edge polynomial is defined as 

i—l N N i-l 

(4.19) e ,(o = - e <w = E = s E d/ ^) - 1 E ^(0- 

fc— fc— z fc— i k—Q 

Let Oft(^) € A°(Q; C ) be expressed as in (4.16), then bh = da^ € A^(Q; Ci) is expressed as 

JV TV JV JV 

(4.20) dE«A(o = E( fl '- a -i) e -® = E fa^M) = E^ e »(^)' 

i=0 i=l i=l i=l 

where <5 is the coboundary operator (4.5), applied to the 0-cochain a( ). It therefore satisfies 
(4.10). For derivations and proofs see [28, 41, 49]. Similar to (4.17), the edge basis-functions are 
constructed such that when integrating (£) over a line segment it gives one for the corresponding 
element and zero for any other line segment, so 

rip ( 1 if i = p 

(4.21) fte,(0 = / e,(0 = 



Equations (4.17) and (4.21) show that indeed we have 1Z1 = Id. The fourth-order Lagrange and 
third-order edge polynomials, corresponding to a Gauss-Lobatto grid with N = 4, are shown in 
Figures 1 and 2. 




Figure 1. Lagrange poly- 
nomials on Gauss-Lobatto- 
Legendre grid. 



Figure 2. Edge polynomials 
on Gauss-Lobatto-Legendre 
grid. 



4.5. Bounded projections and interpolation estimates. The mimetic framework uses La- 
grange, £ HA°(Q), and edge functions, e,(£) £ L 2 A 1 (Q), for the reconstruction, I. Because 
we consider tensor products to construct higher-dimensional interpolation, it is sufficient to show 
that the projection operator is bounded in one dimension. A similar approach was used in [17]. 
Due to the way the edge functions are constructed, there exists a commuting diagram property 
between projection and exterior derivative, 



-> HA° 



-> L 2 A' 



-> 



A?, 



A, 1 , 



-* o, 



which gives, for a £ HA°(Q), the one form d^a = 7r/jda in L 2 A 1 (Q). Lagrange interpolation 
by itself does not guarantee a convergent approximation [25], but it requires a suitably chosen 
set of points, — 1 < £ < £i < • ■ ■ < 6v 1- Here, the Gauss-Lobatto distribution is proposed, 
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because of its superior convergence behaviour. For a e H m A°(Q), the a priori error estimate in 
the 7?A°-norm is given by [18], 

(4.22) \\a - ir h a\\ HA a < Ch l \a\ H i+i A o, I = mm(N, m - 1). 

Equation (4.22) also implies that the projection of zero-forms is stable in the HA°(Q), as is shown 
in the following proposition. 

Proposition 4. [41] For a e HA°(Q) and the projection TTh ■ HA° — > A° , there exists the following 
two stability estimates in H A -norm and H A -semi-norm: 

(4.23) ||7rfcO|| ff A0 < CIHIffAO, 

(4.24) \ir h a\ HA o < C\a\ HA o. 

Now that we have a bounded linear projection of zero forms in one dimension, we can also proof 
boundedness of the projection of one-forms. 

Proposition 5. Let a e HA° and b = da e L? A 1 , then there exists a bounded linear projection 
TT h : L 2 A 1 -> A^, such that 

(4.25) hhb\\ L 2 A i<C\\b\\ L 2 A1 . 

Proof. The proof is based on the result of the previous proposition and the commutation between 
the bounded projection operator and the exterior derivative, Lemma 2, 

Ik/i&II^Ai = |7r/ l da| i 2 A i = |d7T, l a| L 2 A i = \ir h a\ HA a < C\a\ HA a = C||da|| L 2 A i = C||6|| L 2 A i. 

□ 

Propositions 4 and 5 show that the projection 717, is a bounded projection operator, based on 
Lagrange functions and edge functions. As for zero forms using Lagrange interpolation, we can 
also give an estimate for the interpolation error of one forms, interpolated using edge functions. 

Proposition 6. [41] Let a G HA° and b = da G L 2 A 1 , the interpolation error b — 7r^& G L 2 A 1 is 
given by 

(4.26) ||6 - 7r h 6|| L 2 A i < Ch l \b\ H i A i, I = mm(N, m - 1). 

The one dimensional results can be extended to the multidimensional framework by means of 
tensor products. This allows for the interpolation of integral quantities defined on fc-dimcnsional 
cubes. Consider a reference element in R 2 , Q = [— 1,1] 2 . Then the interpolation functions for 
points, lines, surfaces (2D volumes) are given by, 



point : 








line : 




= {em 




surface : 




= ei(OS 


)e J (n). 



The approximation spaces are spanned by combinations of Lagrange and edge basis functions, 
H 1 A (O) DA°(Q; Co) := spaniel 

f N,N f -.N.N 

^A 1 (f7)DAi(Q;C 1 ):=span{(Lg) 1 } xspan{(L«) 2 } 

L 2 A 2 (n) D A 2 (Q; C 2 ) := span , =j • 

For the variables vorticity, velocity and pressure in the VVP formulation of the Stokes problem, 
the /i-convergence rates of the interpolation errors in L 2 A fc -norm become, 

(4.27) 11^-^11^-2 =0(h N+s ), \\u-n h u\\ L2A n-i=0(h N ), \\p - n h p\\ L 2 A n = 0(h N ), 

in case the functions (uj,u,p) are sufficiently smooth, where s = 1 for n = 2 and s = for n > 2. 
The interpolation errors in -ff A fc -norm become, 

(4.28) - n h oo\\ HA n-, = 0(h N ), \\u - Tr h u\\ H A^ - 0(h N ), 
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with N defined as in Section 4.4. 



5. Error estimates 



Next consider the finite dimensional problem: find (oJh,Uh,Ph) € {Wh x V/, x Qh}, given / e 
i 2 A n_1 (ri) and g e £ 2 A™ and boundary conditions in (2.19), for all {r h ,v h ,q h ) E {Wh xV/,x Q^}, 
such that 

(5.1a) a(77,,w/,) - c(r h ,u h ) = h(r ft ), Vr/, e W h , 

(5.1b) e(v h ,u h ) + c(iv h ,Vh) + b( v h,Ph) = f{vh), ^v h eV h: 

(5.1c) b(u h , q h ) = g(qh)- Vqh G Qft- 

The following theorem gives the a priori error estimates of this problem when using the compatible 
spectral discretization method described in the previous section. Corollary 2 showed that we have 
Zw h C Zyy and Zy h C Zy. From this it follows that we have compatible finite dimensional 
subspaces: Wh C W, Vh = dWh © d*Q h C V and Q h = dVh C Q. The derivations of the 
error estimates are based on the methodology of [15]. The proofs are given in the subsequent 
propositions. 

Theorem 2 (Error estimates). Let (uj,u,p) be the solution of the continuous problem given 
in (3.7) or (3.9) and (uJh,Uh,Ph) the solution of the finite dimensional problem in (5.1). The 
continuous problem is well-posed by Propositions 2 and 3 and the finite dimensional problem is 
well-posed by Theorem 1 and Propositions 4 and 5. Furthermore, from Corollary 2 we have that 
for the compatible spectral discretization method, Zw h C Zw and Zy h C Zy . Then the following 
a priori error estimates for the VVP formulation of the Stokes problem hold: 

(5.2) \\u-u>h\\w < + + inf \\w-r h \\ w , 

V a hJ \ IhJ T h £W h 

(5.3) 

\\u~ u h \\ v < ( 1 + — ) ( 1 + ^- J inf ||u - Vh\\v + — ( 1 + — ) ( 1 + — ) inf \\uj-t h \\w, 
V IhJ \ p h Jv h ev h j h \ a h J \ j h J r h ew h 

(5.4) \\p-p h \\ Q <(l+ C f] inf \\ p - qh \\ Q + ^(l+^) U + f) ini \\\u-v h \\ v 

\ PhJ iheQ h p h \ j h J \ p h J v h ev h 

+ (- + 7r)( 1 + -)( 1 + -) ||c-r,|| w . 

\7h PhJ V a hJ V IhJ r h eW h 

Proof. The proof of this Theorem will be given in a series of Propositions 7 to 10. □ 
Proposition 7 (Vorticity error bound). Let ah € Zw h > ^ e error f or vorticity is bounded by 

(5.5) lie*; — LJhWw < ( 1 H — - ) inf lit*; — ahWw- 

V a hJ ° h ez± h 

Proof. Subtract the velocity- vorticity relation in the finite dimensional problem (5.1a) from that 
of the continuous problem (3.9a), we get 

a(r h , u-ujh)- c(jh, u-u h ) = 0, Vr/, e Z Wh C Z w . 

Bound o~h — tOh G Z Wh using inf-sup condition (3.13), we get 

a(T h ,o- h - LU h ) 

a h \\a h - LJh\\w < sup n — n 

r h ez Wh \\Th\\w 

a(rh,o- h - w) + a(T ft ,w - a;/,) 

= sup n — n 

r h ez w , \\Th\\w 



sup 



h 

a(Th,o-h -w) +c(rh,«-Uh) 
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The last term vanishes since r h e Z Wh and Z Wh C Z Wl hence a/Jo^ — < c a|| w — °7i||w- By 

the triangle inequality and the infimum over all ah £ Z^ h we obtain (5.5). □ 

Proposition 8 (Velocity error bound). Let Sh £ Zy h , the error for velocity is bounded by 
(5.6) ||«-« h ||v < fl + -) inf \\u- Sh \\v + -(l + —) inf 



\\0J - <Th\\W- 

Proof. Use the inf-sup condition (3.14) to bound Sh — Uh & Zy, 

c(T h ,S h -U h ) 

~fh\\8h - Uh\\v < sup — n 

c(r /l , s ft - u) + c(r h , u - u h ) 

= sup — n 

r h ez Wh \\ T h\\w 

c(T h ,s h -u) + a(r h ,oj - oj h ) 

= SUp r r 

r h ez Wh \\ T h\\w 
< c c \\u - Sh\\v + c a \\co - U)h\\w- 
By triangle inequality, estimate (5.5) and the infimum over all Sh <G Zy h , we obtain (5.6). □ 
Proposition 9 (Pressure error bound). The error for pressure is bounded by 

\\p-p h \\ Q <{l+ C f] inf\\p- qh \\ Q + C f(l + ^) inf ± \\u-a h \\ v 
V PhJ iheQ h fj h \ j h J Sh ez^ 

(5-7) 

+ (- + 7T 1 + — inf Wu-'hWw- 
\lh PhJ V a hJ <r h ez± h 

Proof. Subtract (5.1b) from (3.9b), we get 

c(u - u h ,v h ) + e(v h ,u - u h ) + b(v h ,p - p h ) = 0, Vv h £ V h . 
So for gj, e Qj, we have 

b(v h ,qh -Ph) = -c(lj - oj h ,v h ) - e(v h ,u- u h ) - b(v h ,p- q h ). 
Use this and the inf-sup condition (3.16) to bound qu — Ph £ Qh, 

b(v h ,q h - Ph) 



Phhh -Ph\\Q < sup 

v h ev h \\Vh\\v 

-c(w - uJh,v h ) - &(v h ,u- u h ) - b(v h ,p- Qh) 

= sup — 

v h ev h \\ v h\\v 

< c c \\uj - uj h \\w + c e \\u - u h \\v + Cb\\p - Q7»||q- 

By triangle inequality, estimates (5.5) and (5.6), and the infimum over all qh £ Qh, we obtain 

(5.7) . □ 

Next we replace the infimums over ah £ ^w h ana - s h € Zy h by best approximation errors. 
Proposition 10. The terms inf rThe z ± \\ UJ ~ cr h\\w o,ndird SheZ ± \\u—Sh\\v are bounded by the best 

approximation estimates iitf The \y h — T~h\\w andmi Vhe v h \\ u ~ v h\\v , using the inf-sup conditions 
(3.16) and (3.14), as 

(5.8) inf \\w-o- h \\ w < + — ) inf \\uj-r h \\ w , 
<Jhez^ h \ jhj r h ew h 

(5.9) inf ||u - SfcHv < (i + lf) inf 

s h ez^ h \ PhJ v h ev h 
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Proof. Take r h <G Wh, then there exists a Kh € Wh such that 

c(K h ,V Zh ) = C(U) - T h ,V Zh ), Vv Zh e 

This is equivalent to 

c^+T/,,^) = c(w,u z J = f{v Zh ), Vv Zh e Zy h , 

which shows that er^ = + t/, e ^viv We can bound ||/t/j||w using the discrete inf-sup condition 
as follows 

c(Kh,u) c(u> - r h ,u) 

7h||«h||w < SUp -jj-jj = Slip jj-j: <c c \\uj-T h \\ v . 

vez Vh \\v\\v vez Vh \\ v \\v 

By triangle inequality and since t/, € Wh was arbitrary, we find (5.8). A similar proof holds for 

(5.9) (see also [15], Proposition 2.5). □ 

Additionally, following section 7.7.6 in [10], we have the following L 2 A k (tt) error estimates for 
the curl of vorticity and divergence of velocity, 

Proposition 11. The errors of the curl of vorticity and divergence of velocity are bounded by 
their best approximation estimates, 

(5.10) ||d(w-w h )|| i2A „-i < inf ||d(w-r h )|| L2A „-i, 

r h ew h 

(5.11) ||d(«-« fc )|| L2A »< mf \\d(u-v h )\\ L 2 An . 

Proof. Choose v = v Zh e Z Vh in (3.9b) and (5.1b) and subtract these. Set v Zh — Ar hl this gives 
the orthogonality relation 

(d(w - u h ), dr h )) n = 0, Vr h e W h - 
Substitute this into the following Cauchy-Schwarz inequality, 

||d(w - ujh)\\ 2 L 2 A „-i = (d(w - uj h ),d(u - T h ))n 

< ||d(w-W/ l )|| L 2 A „-l||d(w-77 l )|| L 2 A n-l ) VT h e W h , 

and (5.10) follows. Next choose qh = dvh € Qh m (3.9c) and (5.1c) and subtract these. Then 

(5.11) follows again from the Cauchy-Schwarz inequality. □ 

Because the projections of respectively to, u, and p, belong to the finite dimensional subspaces 
Wh C W, Vh C V and Qh C Q, the best approximation errors can be bounded using the 
interpolation errors, 

in L \\u- T h\\w < \\u-n h u\\ w , inf \\u-v h \\v < \\u-ir h u\\ v , inf \\p - q h \\ Q < \\p - ir h p\\ Q , 
-r h ew h v h ev h qu^Qu 

and therefore we obtain the following optimal a priori error estimates, 

(5.12) \\u-u;h\\w = 0(h N ), \\u-u h \\v=0(h N ), \\p - p h \\ Q = 0(h N ). 

So the convergence rates for the approximations are equal to those of the interpolations, (4.27), 
(4.28), thus we obtained optimal convergence. The error estimates were obtained independent of 
the chosen types of boundary conditions. 

Remark 2. In contrast to [4, 24] and [10], Table 7.5, where Raviart- Thomas elements were 
used, the proposed compatible method has provably optimal convergence also with standard velocity 
boundary conditions and with non-affine mappings. 
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6. Numerical Results 

Now that the compatible spectral discretization method and its a priori error estimates are 
derived, we perform a series of test problems to show optimal convergence behavior. Purpose of 
the testcases is to show convergence behavior in case of various boundary conditions and in case 
of curvilinear meshes. In all cases we show optimal convergence. 

The first three testcases originate from a recent paper by Arnold et al [4], where suboptimal 
convergence is shown for normal velocity - tangential boundary conditions in vector Poisson and 
Stokes problems, when using Raviart-Thomas elements [48]. Since Raviart-Thomas elements are 
the most popular 77(div, f2) conforming elements, we compare our method to these results. 

6.1. Vector Poisson problems. Figure 3 shows the result of the vector Poisson problem (3.5) 
on fl — [0, l] 2 with coordinates x := (x,y), for a 1-form u G i/A 1 (17), where Y = Y2, i.e. 
with tangential velocity - divergence-free boundary conditions (tr * it = 0, tr * du = 0). The 
corresponding solution is given by 

Vr' = — u(x) dx + u(x) dy 
(6.1) = — (2 sin ttx cos Try) dx + (cos ttx sin Try) dy. 

Both Raviart-Thomas and mimetic spectral clement methods show optimal convergence rates. All 
results of this and the following two problems where obtained on the same quadrilateral mesh of 
2" x 2" subsquares, n = 1, 2, 3, 4, . . . just like the reference solutions from [4]. 




Figure 4 shows again results for the vector Poisson problem for a 1-form, but now in combination 
with normal velocity - tangential velocity boundary conditions (tr u = 0, tr -ku = 0), so Y = Yi. 
The corresponding manufactured solution is 

vP'' = —v(x) dx + m(x) dy 
(6.2) = — (sin ttx sin Try) dx + (sin ttx sin Try) dy. 
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The compatible spectral discretization method again shows optimal convergence, as was expected 
from the above analysis. The Raviart-Thomas elements only show suboptimal convergence in case 
of velocity boundary conditions. This suboptimality was proven in [4], Especially for LOh and 
Auih the current method outperforms the Raviart-Thomas elements, with a difference in rate of 
convergence of |. 




Figure 4. Comparison of the /i-convergence between Raviart-Thomas and 
Mimetic spectral elements for the 2D 1-form Poisson problem with tangential 
velocity - normal velocity boundary conditions. 



6.2. Stokes problems. The same difference in convergence behavior is found for the Stokes 
problem, where T = Ti, i.e. with normal velocity - tangential velocity boundary conditions, see 
Figure 5. Again f2 is the unit square, and the velocity and pressure fields are given by 

w 1 ' = — v(x) dx + u(x) dy 

(6.3) = - (2y 2 (y - l) 2 x(2x ~ l)(x - 1)) dx + (-2x 2 (x - l) 2 y(2y - l)(y - 1)) dy, 

(6.4) p {2) = p(x) dx A dy = ((x ~ \ f + (y - ^) 5 ) dx A dy. 

While for velocity both methods show optimal convergence, for pressure a difference of \ is noticed 
in the rate of convergence and for vorticity and the curl of vorticity again a difference in rate 
of convergence of | is revealed. The error in divergence of velocity is not shown here for the 
Stokes problem, because the method is pointwise divergence-free up to machine precision. Special 
attention to this property is given in [40] . 

We would like to remark is that the results shown in Figure 5 are independent of the kind 
of boundary conditions used. Table 1 shows the results of vorticity for all types of admissible 
boundary conditions. 

The next testcase reveals the optimal convergence in case of higher-order approximation on 
curvilinear quadrilateral meshes for all admissible types of boundary conditions. The manufactured 
solution Stokes problem is given on a curvilinear domain, defined by the mapping (x,y) = <&(£, rj), 
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Figure 5. Comparison of the ^.-convergence between Raviart-Thomas and 
Mimetic spectral element projections for the 2D Stokes problem with normal 
velocity - tangential velocity boundary conditions. 



normal velocity 
tangential velocity 


tangential velocity 
pressure 


vorticity 
normal velocity 


vorticity 
pressure 


convergence 
rate 


1.0280e-04 
1.2445e-05 
1.5424e-06 
1.9238e-07 
2.4035e-08 
Table 1. This U 


1.0109e-04 
1.2410e-05 
1.5426e-06 
1.9247e-07 
2.4042e-08 
ible shows the vortic 


1.0030e-04 
1.2364e-05 
1.5399e-06 
1.9230e-07 
2.4032e-08 
ity error ||w — u)h\ 


1.0035e-04 
1.2375e-05 
1.5416e-06 
1.9255e-07 
2.4065e-08 
^2^0 obtaine 


3.14 
3.05 
3.01 
3.00 
3.00 
d using the 



four types of boundary conditions given in (2.19). The results are obtained on 
an uniform Cartesian mesh with N — 2 and h — | , ^ , ^ , , ^ ■ All four cases 
show third order convergence. 



(6.5a) x(£, v) = h + Ht+To cos(2<) sin^vrr/)) , 

(6.5b) 77) = \ + \ (rj + i sin(27r£) cos^ttt;)) . 

A 6 x 6 element N = 6 mesh is show in Figure 6. Each side of the domain has a different type of 
boundary condition, so T = Ti U T2 U U T4, as shown in the same figure and listed in (2.19). 
The solutions of vorticity lj £ A (51), velocity u £ A 1 (57) and pressure p £ A 2 (57) are given by 

(6.6a) w (0) = §7rsin(§7ra) sin(§7n/), 

(6.6b) — — (cos(§7ra;) sin(|7r?/)) da; + (2 sin(|7ra;) cos(§7ry)) dy, 

(6.6c) p' 2 ^ — (sin(7rx) sin(7ry)) da; A dy. 

They lead to nonzero body force / <E A 1 (57) and mass source g £ A 2 (57). Figure 6 shows the 
convergence of the vorticity uj h £ A° (57; C ), velocity Uh £ A^(57; Ci) and pressure ph £ A 2 t (57; C 2 ). 



22 



JASPER. KREEFT AND MARC GERRITSMA 



The errors for the vorticity and velocity are measured in the i/A fc -norm, i.e. \\lu — u^HhaOj an d 
\\ u ~ w ft II jta 1 j respectively, and the error of the pressure is given in the _L 2 A 2 -norm. In Figure 6 
convergence rates are added which show the optimal /i-convergence behavior of the Stokes problem 
on a curvilinear domain with curvilinear grid and all four types of boundary conditions. 



tangential velocity - pressure 




Figure 6. Upper left figure show the computational domain with a 6 x 6 element 
mesh of N = 6. Furthermore the velocity, vorticity and pressure /i-convergence 
results are shown of Stokes problem (6.6). All variables are tested on grids with 
N = 2, 4, 6 and 8. 

7. Concluding remark 

Optimal approximation of the Stokes problem for all admissible boundary conditions essentially 
hinges on the construction of a conforming discrete Hodge decomposition, A^ = Z% © Z h ' and a 
discrete Poicare inequality, that are based on the bijection of the exterior derivative on the conform- 
ing subspace, d : -2^ ,_L —> B k+1 . Ensuring these properties result in a compatible discretization 
method, and relied on the construction of a bounded projection operator, 7r^ : A k (fl) — > A^(f2; C/c), 
that commutes with the exterior derivative, iv^d = 151Z = diTh- So the compatibility is based on 
the bijection of the coboundary operator, S : Z k,J - —> B k+1 , and the construction of interpolatory 
basis functions. From this it follows that, B^ +1 C B k+1 , Z£ C Z k and Z k ^ C Z k - X . From these 
properties the rest follows. 

For piecewise sufficiently smooth mappings, the optimal conference rates hold on curvilinear 
grids as well, since the pullback operator of the map from a curvilinear domain to the Carte- 
sian frame commutes with the projection operator. Any projection (discretization) with these 
properties will yield similar results as described in this paper. 
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